import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# 1. 数据读取和预处理
data = pd.read_csv('temperature.csv')
print("数据列名:", data.columns.tolist())

# 假设数据包含time和temperature列
time = data['time'].values
temperature = data['temperature'].values

# 创建假设的输入信号
input_signal = np.zeros_like(time)
# 假设在某个时间点有阶跃输入(例如当温度开始上升时)
temp_diff = np.diff(temperature)
step_start_idx = np.argmax(temp_diff > 0.1) + 1 if any(temp_diff > 0.1) else len(time)//10
input_signal[step_start_idx:] = 1  # 假设阶跃输入为1

# 2. 系统辨识
def two_point_method(t, y, u):
    try:
        # 找到输入阶跃变化点
        u_diff = np.diff(u)
        step_idx = np.argmax(u_diff > 0.5) if any(u_diff > 0.5) else len(t)//10
        step_time = t[step_idx]
        step_magnitude = u[step_idx+1] - u[step_idx]
        
        # 相对阶跃时间
        t_rel = t[step_idx:] - step_time
        y_rel = y[step_idx:] - y[step_idx]
        
        # 找到阶跃响应的关键点
        y_final = y_rel[-1]
        if abs(y_final) < 1e-6:
            raise ValueError("系统响应太小")
        
        y_28 = 0.28 * y_final
        y_63 = 0.63 * y_final
        
        # 找到对应的时间点
        t28 = t_rel[np.argmax(y_rel >= y_28)]
        t63 = t_rel[np.argmax(y_rel >= y_63)]
        
        # 计算参数
        tau = 1.5 * (t63 - t28)
        T = t63 - tau
        K = y_final / step_magnitude
        
        return K, tau, T
    except Exception as e:
        print(f"两点法失败: {str(e)}")
        return None

# 尝试系统辨识
ident_result = two_point_method(time, temperature, input_signal)
if ident_result:
    K, tau, T = ident_result
    print(f"辨识结果: K={K:.2f}, tau={tau:.2f}s, T={T:.2f}s")
else:
    # 使用图片中提供的默认参数
    K = 0.99
    tau = 2895
    T = 0
    print("使用默认参数:", f"K={K}, tau={tau}, T={T}")

# 3. PID控制器类
class PIDController:
    def __init__(self, Kp, Ki, Kd, setpoint, output_limits=(0, 100)):
        self.Kp = Kp
        self.Ki = Ki
        self.Kd = Kd
        self.setpoint = setpoint
        self.output_limits = output_limits
        self.reset()
    
    def reset(self):
        self.last_error = 0
        self.integral = 0
        self.last_time = None
    
    def __call__(self, measurement, dt):
        error = self.setpoint - measurement
        
        # 比例项
        P = self.Kp * error
        
        # 积分项
        self.integral += error * dt
        I = self.Ki * self.integral
        
        # 微分项
        if self.last_time is not None and dt > 0:
            derivative = (error - self.last_error) / dt
            D = self.Kd * derivative
        else:
            D = 0
        
        # 总和
        output = P + I + D
        
        # 限幅
        output = np.clip(output, *self.output_limits)
        
        # 保存状态
        self.last_error = error
        self.last_time = (self.last_time or 0) + dt
        
        return output

# 4. 闭环系统仿真和优化
def simulate_system(pid_params, model_params, setpoint, initial_temp, sim_time=10000):
    Kp, Ki, Kd = pid_params
    K, tau, T = model_params
    
    pid = PIDController(Kp, Ki, Kd, setpoint)
    t_sim = np.linspace(0, sim_time, 1000)
    dt = t_sim[1] - t_sim[0]
    y = np.zeros_like(t_sim)
    u = np.zeros_like(t_sim)
    y[0] = initial_temp
    
    for i in range(1, len(t_sim)):
        u[i] = pid(y[i-1], dt)
        # 使用模型模拟系统响应
        if t_sim[i] >= T:
            dy_dt = (K*u[i] - y[i-1]) / tau
            y[i] = y[i-1] + dy_dt * dt
        else:
            y[i] = y[i-1]
    
    return t_sim, y, u

def evaluate_performance(t, y, setpoint):
    error = setpoint - y
    rise_time_idx = np.argmax(y >= 0.9 * setpoint) if any(y >= 0.9 * setpoint) else len(t)-1
    rise_time = t[rise_time_idx]
    
    settling_mask = np.abs(error) < 0.02 * setpoint
    settling_time_idx = np.argmax(settling_mask) if any(settling_mask) else len(t)-1
    settling_time = t[settling_time_idx]
    
    overshoot = max(0, np.max(y) - setpoint)
    steady_state_error = np.mean(error[-100:])  # 最后100个点的平均误差
    
    return {
        'rise_time': rise_time,
        'settling_time': settling_time,
        'overshoot': overshoot,
        'steady_state_error': steady_state_error
    }

# 5. 主程序
setpoint = 35.0  # 设定值35℃
initial_temp = 16.8  # 初始温度

# 运行仿真并绘制结果
t_sim, y_sim, u_sim = simulate_system([1.0, 0.001, 10.0], (K, tau, T), setpoint, initial_temp)

# 评估性能
performance = evaluate_performance(t_sim, y_sim, setpoint)
print("\n初始PID参数性能:")
for k, v in performance.items():
    print(f"{k}: {v:.2f}")

# 绘制结果
plt.figure(figsize=(12, 6))
plt.plot(t_sim, y_sim, label='温度响应')
plt.axhline(setpoint, color='r', linestyle='--', label='设定值')
plt.title('温度控制系统响应')
plt.xlabel('时间 (s)')
plt.ylabel('温度 (℃)')
plt.legend()
plt.grid(True)
plt.show()
